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Abstract 

We  develop  a phase-field  model  for  the  solidification  of  a pure  material  that  includes 
convection  in  the  liquid  phase.  The  model  permits  the  interface  to  have  an  anisotropic 
surface  energy,  and  allows  a quasi-incompressible  thermodynamic  description  in  which 
the  densities  in  the  solid  and  liquid  phases  may  each  be  uniform.  The  solid  phase  is 
modeled  as  an  extremely  viscous  liquid,  and  the  formalism  of  irreversible  thermodynamics 
is  employed  to  derive  the  governing  equations.  We  investigate  the  behavior  of  our  model 
in  two  important  simple  situations  corresponding  to  the  solidification  of  a planar  interface 
at  constant  velocity:  density  change  flow  and  a shear  flow.  In  the  former  case  we  obtain  a 
non-equilibrium  form  of  the  Clausius-Clapeyron  equation  and  investigate  its  behavior  by 
both  a direct  numerical  integration  of  the  governing  equations,  and  an  asymptotic  analysis 
corresponding  to  a small  density  difference  between  the  two  phases.  In  the  case  of  a parallel 
shear  flow  we  are  able  to  obtain  an  exact  solution  which  allows  us  to  investigate  its  behavior 
in  the  sharp  interface  limit,  and  for  large  values  of  the  viscosity  ratio. 
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1 Introduction 


Diffuse  theories  of  interfaces  separating  two  bulk  phases  were  developed  in  the  19th  Century  by 
Poisson  [1],  Gibbs  [2],  Maxwell  [3],  Rayleigh  [4],  and  van  der  Waals  [5].  Previously,  interfaces 
had  been  modeled  by  Young,  Laplace  and  Gauss  as  surfaces.  In  the  latter  formulation  the 
interface  is  regarded  as  a singular  surface  on  which  associated  physical  mechanisms  are  localized 
and  represented  as  boundary  conditions  to  be  applied  at  the  surface;  e.g.,  the  notion  of  surface 
energy  as  an  energy  per  unit  area  of  the  interfacial  surface.  This  description  of  a phase  boundary 
is  sometimes  referred  to  as  a ‘sharp  interface’  model  and  results  in  a so-called  ‘free  boundary 
problem’.  In  contrast,  diffuse  interface  theories  recognize  that,  in  reality,  the  interface  has  a 
finite  thickness  (albeit  small  compared  with  typical  macroscopic  length  scales)  in  which  physical 
quantities,  such  as  density  or  composition,  vary  between  their  values  in  the  adjacent  bulk  phases 
(see,  e.g.,  [6,7]).  Quantities  that  in  the  sharp  interface  formulation  are  regarded  as  localized  to 
the  interfacial  surface  are,  in  the  diffuse  interface  setting,  identified  as  being  distributed  within 
the  interfacial  region.  For  example,  the  surface  energy  of  an  isothermal  interface  is  derived  from 
the  elevated  Helmholtz  free  energy  density  throughout  the  whole  interfacial  region. 

Diffuse  interface  models  may  be  based  on  an  extended  thermodynamics  involving  gradients 
of  the  thermodynamic  variables  to  account  for  nonlocal  effects.  Originally  such  theories  were 
formulated  to  investigate  liquids  near  their  critical  point  and  have  subsequently  been  refined 
and  developed  to  account  for  a wide  range  of  physical  situations,  such  as  liquid  crystals  [8], 
superconductivity  [9],  spinodal  decomposition  [10,11]  and  ordering  transitions  in  alloys  [12-14]. 
Rowlinson  and  Widom  [15]  provide  a thorough  account  of  their  historical  development. 

The  phase-field  model  of  the  first-order  phase  transition  associated  with  the  solidification  of 
a pure  material  was  first  proposed  by  Langer  [16, 17]  and  subsequently  developed  by  a number 
of  researchers  [18-23].  Phase-field  models  provide  an  example  of  a diffuse  interface  model  in 
which  an  order  parameter,  <j>,  is  postulated  whose  value  indicates  the  phase  of  the  system  at  a 
particular  point  in  space  and  time  (in  this  paper  <j>  = 1 and  <f>  = 0 denote  the  solid  and  liquid 
phases,  respectively).  Langer  represented  the  free  energy  of  a single-component  system  by  a 
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gradient  energy  functional  of  the  form 


T = Jv{lt2\V<f>\2  + f(4>,T)}dV, 


(1) 


where  e is  the  gradient  energy  coefficient  and  T is  the  temperature.  The  free  energy  density, 
/(0,  T),  has  a double- well  structure  with  respect  to  0 in  which  the  two  local  minima  correspond 
to  the  solid  and  liquid  phases.  Langer  proposed  the  following  governing  equations  for  the  phase 
field  and  temperature: 


,,<90  6T 

M—  = — — - 

dt  5(f) 


(2) 

(3) 


where  1/M  is  a positive  constant  termed  the  mobility,  c is  the  heat  capacity,  k is  the  ther- 
mal conductivity  and  L is  the  latent  heat  per  unit  volume  of  the  material.  This  phase-field 
formulation  replaces  the  free-boundary  problem  associated  with  the  sharp  interface  model  of 
an  interface  by  a coupled  pair  of  nonlinear  reaction  diffusion  equations.  The  location  of  the 
interface  is  represented  by  the  level  set  0 = 1/2. 

The  original  derivation  of  the  phase-field  equations  was  justified  by  requiring  the  free  energy 
of  the  system  to  decrease  monotonically  in  time.  The  equation  for  the  temperature  (3)  wras 
based  on  a modification  of  the  heat  equation  to  allow  a source  term  that  accounts  for  latent  heat 
production  at  a moving  interface.  Subsequently,  Penrose  and  Fife  [23]  and  others  [24-27]  applied 
the  arguments  of  irreversible  thermodynamics  to  the  derivation  of  the  phase-field  equations, 
establishing  that  they  are  consistent  with  non-negative  local  entropy  production. 

The  phase-field  equations  have  been  extended  to  allow  for  anisotropic  surface  energy  by  a 
number  of  authors  [20,28,29].  In  particular,  Kobavashi  [28]  proposed  that  the  gradient  energy 
coefficient,  e,  may  be  regarded  as  a function  of  V0  in  order  to  model  surface  energy  anisotropy, 
and  Taylor  and  Cahn  [30,31]  suggested  that  the  square  gradient  term  in  the  free  energy  functional 
be  replaced  by  [T(V0)]2,  where  T(V0)  is  a homogeneous  degree  one  function  of  its  argument. 
Wheeler  and  McFadden  [32]  showed  how  this  formulation  could  be  used  to  define  a generalized 
form  of  the  ^-vector  [33,34],  which  provides  an  elegant  description  of  surface  energy  anisotropy. 
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If  we  denote  p = V</>,  then  the  f-vector  is  a homogeneous  function  of  p of  degree  zero  with 
components  given  by 

tj(p)  = £-r  (p);  (4) 

dpj 

£ satisfies  the  fundamental  relation  [33] 

dT(p)=£-dp.  (5) 

For  an  isotropic  surface  energy,  this  gives  £ = V</>/|V0|. 

The  relation  between  the  phase-field  equations  and  the  free-boundary  formulation  can  be 
established  by  examining  the  so-called  sharp  interface  limit  of  the  phase-field  equations.  The 
phase-field  equation  has  solutions  in  which  fronts  form  with  a width  proportional  to  e wherein 
<j>  varies  between  zero  and  one;  the  fronts  represent  interfaces.  By  investigating  solutions  of  the 
phase-field  equations  in  the  limit  e — > 0,  Caginalp  [22]  showed  that  the  governing  equations  of  the 
phase-field  model  converge  to  a sharp  interface  model  with  associated  boundary  conditions,  i.e., 
a free-boundary  problem.  By  considering  different  distinguished  limits  Caginalp  showed  that  a 
variety  of  different  free-boundary  problems  emerge  in  the  sharp  interface  limit,  most  of  which  are 
different  forms  of  the  classical  Stefan  problem,  as  well  as  the  Hele-Shaw  problem.  Karma  and 
Rappel  [35]  also  examined  different  distinguished  limits  that  allow  more  flexible  interpretations 
of  interface  kinetics  to  be  used  in  numerical  calculations.  For  the  anisotropic  phase-field  model 
McFadden  et  al.  [36]  obtained,  in  the  sharp  interface  limit,  the  anisotropic  version  of  the  Gibbs- 
Thomson  equation  in  two  dimensions.  Subsequently,  Wheeler  and  McFadden  [32]  extended  this 
analysis  to  three  dimensions  by  employing  the  generalized  f-vector. 

One  of  the  main  advantages  of  the  phase-field  formulation  over  the  free-boundary  problem 
lies  in  its  potential  to  compute  realistic  and  complex  interface  shapes  associated  with  dendritic 
growth.  Early  calculations  [37]  were  restricted  by  the  available  computing  power  to  simple 
interface  geometries.  Kobayashi  [28]  identified  the  importance  of  surface  energy  anisotropy 
for  the  computation  of  dendritic  growth  and  exhibited  computational  results  corresponding  to 
both  two  and  three-dimensional  dendrites.  Subsequently,  numerous  workers  [29,  38-41]  have 
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refined  the  computational  technique,  provided  numerical  solutions  with  improved  accuracy,  and 
addressed  issues  related  to  the  viability  of  the  computational  approach. 

The  phase  field  model  described  above  does  not  include  coupling  to  the  momentum  equation 
and  viscous  stress  tensor  of  classical  hydrodynamics.  The  aim  of  the  present  paper  is  to  formulate 
a phase-field  model  that  is  fully  coupled  to  the  equations  of  hydrodynamics  in  a self-consistent 
manner.  We  shall  find  that  this  requires  the  inclusion  of  an  additional  stress  tensor  associated 
with  the  diffuse  interface.  Wheeler  and  McFadden  [42]  showed,  using  Noether’s  theorem  (see, 
e.g.,  [43])  that  associated  with  the  steady  phase-field  equation  (2)  there  exists  a second  rank 
tensor 


-e2|V0|2  + / 


I — e2W(t)  0 V0, 


(6) 


where  0 denotes  the  outer  product  (see  Appendix)  and  I is  the  unit  tensor,  which  satisfies  the 
conservation  lawr  V • E = 0.  E represents  the  part  of  the  stress  tensor  that  results  in  capillary 
forces  within  the  interfacial  region.  This  conservation  law  has  also  be  derived  by  Fried  and 
Gurtin  [44,45]  using  an  alternative  approach  employing  configurational  forces.  The  tensor  S is 
the  counterpart  of  the  capillary  tensor  [46]  that  acts  as  the  reversible  part  of  the  stress  tensor 
in  the  theory  of  fluids  near  a critical  point;  the  irreversible  part  is  provided  by  the  standard 
viscous  stress  term  of  a Newtonian  liquid. 

In  the  original  theory  for  a critical  fluid,  the  density,  which  satisfies  the  continuity  equation, 
is  treated  as  the  order  parameter  and  appears  instead  of  4>  in  the  capillary  tensor  analogous  to 
(6).  The  momentum  equation,  modified  to  include  the  divergence  of  the  capillary  tensor,  governs 
the  flow  while  an  equation  of  state  relates  the  pressure,  temperature,  and  density.  This  original 
theory  for  a critical  fluid  has  been  extended  to  investigate  a range  of  hydrodynamic  phenomena 
including  capillary  waves,  moving  contact  lines,  droplets  and  nucleation  [47].  In  the  context  of 
a binary  fluid,  the  composition  may  play  the  role  of  a conserved  order  parameter  that  satisfies 
a Cahn-Hilliard  equation  [10].  A variety  of  situations  have  been  studied  ranging  from  spinodal 
decomposition  to  thermocapillary  flow;  the  review  by  Anderson  et  al.  [47]  and  references  therein 
provide  further  details.  An  early  attempt  to  include  fluid  motion  within  a phase-field  model 


5 


is  due  to  Caginalp  and  Jones  [48,49].  They  appended  the  inviscid  momentum  equation  and 
the  continuity  equation  to  the  phase-field  model,  but  did  not  address  the  issues  of  momentum 
balance  in  the  solid  and  capillary  contributions  to  the  stress  tensor.  Diepers  et  al.  [50]  have 
employed  the  methodology  of  two-phase  fluid  flow,  where  (j>  is  interpreted  as  a solid  fraction. 
Their  model  is  used  to  study  coarsening  in  a binary  solid/liquid  mixture  with  and  without  flow. 
Tonhardt  and  Amberg  have  also  performed  two-dimensional  numerical  studies  using  adaptive 
finite  elements  to  study  the  effects  of  a shear  flow  on  dendritic  growth  morphology  [51,52]. 

In  this  paper  we  bring  together  several  ideas  to  develop  a phase-field  model  which  allows  for 
convection  in  the  liquid  phase.  Our  model  has  two  notable  aspects:  first,  we  represent  both  the 
solid  and  liquid  phases  as  Newtonian  fluids  in  which  the  viscosity  of  the  putative  solid  phase  is 
specified  to  be  much  larger  than  that  of  the  liquid  phase.  Second,  the  interface  is  ascribed  an 
anisotropic  surface  energy,  which  is  non-standard  for  a model  which  treats  the  two  phases  as 
Newtonian  fluids.  These  unconventional  features  are  in  keeping  with  our  intention  to  model  a 
solid-liquid  system.  In  order  to  obtain  the  desired  viscosity  variation  between  the  phases,  the 
viscosity  is  assumed  to  depend  on  the  phase  field,  (p.  The  anisotropic  surface  energy  is  achieved 
by  employing  the  generalized  £- vector  formalism  [42].  Unlike  previous  diffuse  interface  models, 
which  incorporate  fluid  motion  coupled  to  a conserved  order  parameter  description  [47],  we  adopt 
a nonconserved  order  parameter,  </>,  in  line  with  our  aim  of  directly  extending  conventional  phase 
field  models  of  solidification  to  account  for  convection.  This  has  the  advantage  that  we  may 
treat  quasi-incompressible  systems  [53]  in  which  the  density  of  the  solid  and  liquid  bulk  phases 
are  each  spatially  uniform  by  allowing  the  density,  p,  to  be  a prescribed  function  of  (p. 

With  these  assumptions  we  develop  the  irreversible  thermodynamics  of  the  model  from  gra- 
dient functionals  for  both  the  entropy  and  internal  energy.  We  identify  governing  equations  that 
are  consistent  with  the  first  and  second  laws  of  thermodynamics.  The  quasi-incompressibility 
assumption,  which  allows  the  density  to  depend  solely  on  <p  and  not  the  pressure,  p,  restricts 
the  form  of  the  thermodynamic  potentials  that  may  be  employed  [53].  The  model  comprises  the 
compressible  Navier-Stokes  equations  with  a modified  stress  tensor  which  includes  additional 
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terms  related  to  S as  given  in  equation  (6),  an  energy  equation  and  a phase- field  equation 
involving  a material  time  derivative  of  0. 

In  order  to  investigate  and  validate  this  new  model  we  examine  it  in  some  simple  situations. 
First,  we  consider  an  isothermal  planar  interface  at  equilibrium.  We  show  that  the  model 
recovers  the  classical  equilibrium  state  in  which  the  chemical  potential  and  pressure  are  equal 
in  each  bulk  phase  and  the  changes  in  temperature  and  pressure  for  coexistence  are  related 
to  the  density  and  entropy  mismatch  of  the  two  phases  by  the  Clausius-Clapeyron  equation. 
We  then  investigate  the  non-equilibrium  situation  in  which  the  planar  interface  advances  with 
constant  velocity  and  an  advective  flow  is  induced  in  the  liquid  due  to  the  density  mismatch.  We 
derive  a non-equilibrium  form  of  the  Clausius-Clapeyron  equation,  obtain  numerical  solutions 
and  conduct  an  asymptotic  analysis  for  small  density  mismatch.  Finally,  we  investigate  the 
situation  in  which  a planar  interface  advances  with  constant  speed  into  the  melt  with  a parallel 
shear  flow  ahead  of  it.  We  are  able  to  obtain  closed  form  solutions  and  investigate  them  in  both 
the  sharp-interface  limit  as  well  as  the  limit  in  which  the  viscosity  of  the  solid  is  much  greater 
than  that  in  the  liquid. 

2 The  Model 

We  consider  a non-isothermal  system  consisting  of  a pure  material  that  may  exist  in  two  distinct 
phases.  We  follow  the  standard  phase-field  methodology  and  introduce  a phase-field  variable, 
0(x,  £),  whose  value  indicates  the  thermodynamic  phase  of  the  system  as  a function  of  time,  t, 
and  position,  x.  Both  phases  are  treated  as  fluids,  although  in  the  applications  we  will  assume 
that  one  phase  has  a much  larger  viscosity  and  interpret  it  as  an  approximation  to  a solid  phase. 
In  many  solidification  applications,  a fluid  model  is  used  for  the  thermodynamic  description 
of  the  the  solid  phase,  in  that  the  elastic  properties  of  the  solid  are  ignored.  We  will  also 
consider  that  the  phase  transition  is  first-order  and  has  an  anisotropic  surface  energy,  which  is 
unconventional  for  a fluid-fluid  system,  but  is  consistent  with  our  intention  to  model  a solid- 
liquid  system.  We  adopt  the  convention  that  0 = 0 denotes  the  liquid  phase  and  0 = 1 denotes 
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the  solid  phase.  A solid-liquid  interface  is  represented  by  a thin  layer  in  which  the  phase  field 
varies  rapidly  between  zero  and  unity.  The  governing  equations  are  derived  by  following  the 
formalism  of  irreversible  thermodynamics  [54],  as  originally  applied  to  the  phase-field  equations 
by  Penrose  & Fife  and  others  [23-25].  The  steady-state  versions  of  the  equations  can  also  be 
obtained  by  appealing  to  variational  arguments  as  in  Ref.  [42].  Derivations  based  on  mechanical 
microforce  balance  laws,  as  developed  by  Gurtin  et  al.  [55,56],  are  also  possible. 


2.1  Nonequilibrium  Equations 


We  assume  that  the  total  entropy,  in  a material  volume,  Cl(t),  of  the  system  is  given  by 


<S  = 


= / 

JQ{t) 


ps 


y|r2(V0) 


dV, 


(7) 


where  p is  the  density  and  s is  the  entropy  per  unit  mass.  The  first  term  in  the  integrand,  ps,  is 
the  classical  entropy  density  (per  unit  volume)  and  the  second  is  a nonclassical  term  associated 
with  spatial  gradients  of  the  phase  field.  Here  the  gradient  entropy  coefficient  es  is  assumed 
to  be  a constant  for  simplicity,  and  F is  a homogeneous  function  of  degree  unity.  As  we  show 
below,  the  function  T allows  for  a general  anisotropic  surface  energy  of  the  solid-liquid  interface. 
An  isotropic  surface  energy  results  from  the  choice  T(V0)  = |V</>|. 

The  total  mass,  A4,  linear  momentum,  V,  and  internal  energy,  E,  associated  with  the  material 
volume  are  assumed  to  have  the  form 


M 

V 

E 


dV, 


(8) 

(9) 

(10) 


respectively1.  Here  u is  the  velocity,  e is  the  internal  energy  density  (per  unit  mass)  and  ce  is 

1For  simplicity  we  omit  a possible  gradient  term  in  the  density  functional  M ; its  inclusion  modifies  the 
continuity  equation  and  replaces  the  gradient  Helmholtz  free  energy  coefficient  ep  in  Eqns.  (29)  and  (35)  by 
a gradient  coefficient  ep  corresponding  to  a thermodynamic,  or  Kramers,  potential  function.  This  is  briefly 

indicated  in  section  2.2  below. 
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the  gradient  energy  coefficient,  which  is  assumed  to  be  constant.  The  thermodynamic  relations 

de  = T ds  + dp  + d$,  (11  j 

p 2 d(p 

e = Ts  - p/ p + p,  (12) 

are  assumed  to  apply  locally,  where  p is  the  thermodynamic  pressure  and  p is  the  chemical 
potential  (or  Gibbs  free  energy  per  unit  mass). 

The  physical  balance  laws  for  mass,  linear  momentum,  and  internal  energy  are  given  by 


dM 

dt 

- 0, 

(13) 

de 

dt  . 

dV 

dt 

= / fi  • m dA , 

J 6Q(t) 

(14) 

/ qE  • n dA 
J6n(t) 

II 

oT"  -s 

jo 

3 

Si 

s. 

> 

(15) 

respectively,  where  h is  the  outward  unit  normal  to  <5Q(t) , m is  the  stress  tensor,  and  q%  is 
the  internal  energy  flux.  The  momentum  balance  (14)  requires  that  the  rate  of  change  of  the 
total  momentum  of  the  material  volume  results  from  forces  acting  on  its  boundary  5Q(t)  (for 
simplicity  we  neglect  body  forces  such  as  gravity;  their  inclusion  is  straightforward).  The  energy 
balance  (15)  equates  the  rate  of  change  of  the  total  internal  energy  of  f2(t)  plus  the  energy  flux 
through  its  boundary  to  the  rate  of  w^ork  of  the  forces  at  its  boundary. 

In  addition,  the  entropy  balance  takes  the  form 

— +/  qs  • h dA  = f sprod  dV,  (16) 

dt  Jsn(t)  Jn{t) 

wrhere  qs  is  the  entropy  flux  and  sprod  is  the  local  rate  of  entropy  production.  The  second  law 
of  thermodynamics  is  then  expressed  by  the  requirement  that  sprod  is  non-negative. 

To  proceed  we  recast  the  conservation  laws  ( 13)— (16)  as  differential  equations.  These  are 
used  to  express  the  local  entropy  production  in  terms  of  the  fluxes  m,  and  qs , as  well 
as  D<j)/Dt.  We  then  identify  forms  for  these  quantities  vrhich  ensure  that  the  local  entropy 
production  is  non-negative.  The  fluxes  that  result  from  this  procedure  involve  both  classical 
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contributions  and  non-classical  contributions  that  depend  on  V</>.  In  addition,  we  obtain  an 
evolution  equation  for  the  phase  field. 

Applying  the  Reynolds  transport  theorem  [57]  to  the  mass  balance  law  (13)  gives  the  conti- 
nuity equation  in  its  conventional  form 

^ + PV-u  = 0,  (17) 


where  the  material  derivative  of  p is  denoted  by  Dp/Dt  = dp/dt  + (u  • V)p.  Similarly,  the  linear 
momentum  equation  takes  the  form 


Du 


P 


Dt 


V • m. 


(18) 


where,  if  rrijk  denotes  the  components  of  m,  then  V • m has  components  drrijk/dxj. 

The  energy  equation  is  more  complicated,  and  we  discuss  it  in  more  detail.  From  (15)  it 
follows  that 


/ 

Jn{t) 


De  _ Du  _ 

PlDt  + PU  ' ~Dt  ~ V ' ’ U'  + V * QE 


dV  + if  h|r 2(v<p)dv  = o. 

at  Jn(t)  2 


(19) 


In  Appendix  A we  show  that 


-! 

dt  Jn(t) 


l-r2(V4>)dv  = 


(20) 


where 

Qg  = v • (rf|T  - • (r|)  - rvu : £ ® v^»  + ir2v  • u.  (21) 

Here  we  have  introduced  the  Cahn-Hoffman  ^-vector  [33]  for  a diffuse  interface  [42]  as  given  by 
Eqn.  (4).  If  we  set  p = V0,  it  has  components  £j  = dT (p)/dpj]  for  an  isotropic  surface  energy, 
this  gives  f = V<j>/\V<j)\.  Hence  we  deduce  that 

Hg  J~^) 

p-—  + pu  • — - - V • (m  ■ u)  + V • qE  + e% QG  = 0.  (22) 


We  now  employ  the  conservation  of  linear  momentum  (18)  to  rewrite  this  equation  as 


De 

p — + u • (V  • m)  - V • (m  • u)  + V • qE  + £2eQg  = 0, 


(23) 
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which  on  using  the  identity 


V • (m  • u)  — (V  • m)  • u + m : Vu, 


(24) 


simplifies  to  yield  the  energy  equation  as 


De 

~Dt 


p—  + V • qE  = m : W - e2EQG. 


(25) 


Here  m : Vu  = nrikjduj/ dxk  (with  summation  over  repeated  indices  implied).  In  an  analogous 
way,  the  entropy  balance  (16)  leads  to  the  result 


Ds 

Dt 


Pttt  + V • qs  — sprod  + e| QG. 


(26) 


The  thermodynamic  relation  between  Ds/Dt  and  De/Dt  follows  from  Eqn.  (11)  and  is  given 


by 


De  _ Ds  p Dp  de  D6 
Dt  Dt  p2  Dt  dcj)  Dt 


(27) 


The  continuity  equation  (17)  and  Eqns.  (21),  (25),  and  (27)  can  be  used  to  express  the  entropy 
production  given  by  Eqn.  (26)  as 


\prod  _2_ 

T 


rrt  + e|T£  0 V0  + 


6 Ft  2 


"-fr 


+ V • ( qs  - ) + ( & + 4r(^)  ■ v ( ^ ) , 


& _ eF  TC0^ 

T T ^ Dt 


Dt  J 


(28) 


where  e2F  = e2E  + Te|. 


We  now  make  the  following  choices  for  the  constitutive  equations  for  the  fluxes  and  D<j>/Dt 
which  ensure  that  sprod  is  positive 


M 


m = 

D(j) 

~Dt 

Qe  = 

Qs  = 


-P  + fr- 


I - cjTf  0 V0  + r, 


*(?)-**& 


(29) 

(30) 

(31) 

(32) 
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Here  r is  the  viscous  stress  tensor,  which  for  a Newtonian  fluid  is  given  by  t — p(Vu  + S7ur)  + 
A(V  • u)I , where  p and  A are  coefficients  of  viscosity,  I is  the  unit  tensor,  and  M is  a positive 
mobility  coefficient  which  we  take  to  be  constant.  A constant  value  for  the  thermal  conductivity 
k corresponds  to  the  choice  k = T2  k.  The  resulting  equations  of  motion  assume  the  form 


P 


M 


Dp 

Dt 

Du 

'Dt 

D<\> 

~Dt 

De 

Dt 


= -pV  • u, 


V • 


-p  + ^4  r2 ) / - <g>  v<t>  + r 


= 4v-(r  a-P 


de 


= V • [hVT]  + e|V  • (T|) 


Dcf) 

~Dt 


+ 


-P  + jTe |T2  ) I - Te\Yi  ® + r 


: Vu. 


(33) 

(34) 

(35) 

(36) 


2.2  Equilibrium  Equations 


The  equilibrium  form  of  the  above  governing  equations  with  u = 0 admits  an  isothermal  solution 
that  also  satisfies 


0 = V- 


i-p+j  r2)/-4rf® 

o = 4v.(rj)-p|. 


(37) 

(38) 


By  applying  the  divergence  operator  in  (37)  and  using  the  fundamental  relation  (5),  Eqn.  (37) 
can  be  reduced  to  the  form 


o = Vp  + e2FV  ■ (if)  V<f>. 


(39) 


These  equations  also  result  directly  from  a variational  formulation,  which  may  be  stated  equiva- 
lently as  either  an  entropy  maximization  or  an  energy  minimization;  we  briefly  sketch  the  latter 
argument.  We  temporarily  include  a gradient  coefficient  €m  in  the  density  functional  (8)  to  illus- 
trate the  appearance  of  the  thermodynamic,  or  Kramers,  potential  function  in  the  formulation; 
we  subsequently  set  cm  = 0 to  recover  the  simpler  model  that  we  consider  in  the  remainder  of 
the  paper. 
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We  consider  a convection-free  equilibrium  state,  and  set 


0 = S{£  - XsS  - A mM)  = 6 J 


pe  — X sps  — XmP  + 


(ei  + Xse2s  - Am^)  2 


r2(v<« 


dV,  (40) 


where  A 5 and  Am  are  Lagrange  multipliers  that  are  introduced  to  account  for  the  constraints  of 
constant  total  entropy  and  mass.  Taking  variations  with  respect  to  5s,  5p,  and  5</>  gives 


0 = pes-  Xsp , (41) 

0 = e + pep  — XsS  — X m,  (42) 

0 = pe^  — (eE  + A 565  — XMe2M)V  • (T<f),  (43) 

respectively.  Eqn.  (41)  implies  that  the  temperature  T = es  is  uniform  and  equal  to  the  Lagrange 
multiplier  A$.  Using  the  relation  ep  = p/p2,  Eqn.  (42)  gives  that  Am  = e + p/p  — Ts  = p is 
also  constant.  In  Eqn.  (43)  it  follows  that  e2E  4-  Xse2s  — X m^m  = + Te2s  — pe2M  = e2K,  where 


eE  is  a gradient  coefficient  corresponding  to  the  thermodynamic,  or  Kramers,  potential  function 
pe  — pTs  — pp.  Setting  cm  — 0 replaces  eE  by  the  Helmholtz  free  energy  coefficient  eE  defined 
by  e2F  = e\  + A 565,  so  that  Eqn.  (43)  becomes  identical  to  Eqn.  (38).  We  henceforth  set  cm  = 0, 
so  that  eE  reduces  to  ep. 

The  equivalence  of  Eqn.  (43)  and  Eqn.  (39)  results  from  the  Gibbs-Duhem  equation,  which 
for  constant  values  of  T and  p implies  that 

-Vp  = pe^V^  (44) 

Alternatively,  an  application  of  Noether’s  theorem  [43]  to  the  translation  invariant  functional 
C — — Q+  \t2FY2  in  Eqn.  (40),  where  — p = pe  — pTs  — pp,  leads  to  a divergence-free  stress  tensor 
[42] 

<«> 

which  is  identical  to  the  equilibrium  form  for  m in  Eqn.  (29)  with  r = 0 that  appears  in 
Eqn.  (37). 
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Equilibrium  is  therefore  characterized  by  uniform  temperature  and  chemical  potential.  In 
contrast,  there  are  pressure  gradients  in  the  interfacial  regions  where  V</>  / 0,  although  the 
far-held  pressures  in  the  bulk  phases  are  equal. 


2.3  Quasi-Incompressible  Thermodynamics 


In  order  to  study  situations  in  which  the  density  in  each  phase  is  uniform,  it  is  convenient  to 
adopt  a thermodynamic  formation  which  does  not  employ  the  density  as  an  independent  variable, 
as  in  the  model  of  quasi-incompressible  how  considered  by  Lowengrub  and  Truskinovsky  [53] 
and  also  in  the  work  of  Rooney  et  al.  [58]  on  modeling  thermal  expansion  in  a Newtonian  huid. 
We  therefore  choose  the  pressure  and  temperature  as  independent  variables,  and  work  with  a 
Gibbs  free  energy  per  unit  mass,  g(T,p:  </>),  which,  aside  from  its  argument,  is  formally  identical 
to  the  chemical  potential  p(p,T,4>)  appearing  in  Eqn.  (12).  The  internal  energy  per  unit  mass 
may  then  be  written  in  the  form 


e = g{T,p,<t>)  + Ts{T,p,  <j>) 


P 


p{T,p , 4>Y 


(46) 


where  we  note  the  identities 
s(T,p,  <j>)  = -■ 


dg 

<^>  1 

1 

r-H 

de 

. d9 

dT 

p{T,p,<t> ) dp 

T,cp 

dcj) 

S,p  T 

(47) 


V,T 


where  the  variables  which  are  held  constant  in  forming  the  various  derivatives  are  indicated 
explicitly.  We  prescribe  the  density  as  a function  of  the  phase-held  variable  <j>  alone, 


p(4>)  = psr(4>)  + pL[  1 - r(<t>)\ , (48) 

where  the  solid  and  liquid  densities  ps  and  pL  of  the  bulk  phases  are  constants,  and  r(</>)  is  a 
smooth  monotonic  function  that  has  r(0)  = 0 and  r(l)  = 1;  we  will  take  r(</>)  = <j>2{ 3 — 2 </>). 
Then,  from  Eqn.  (47)  we  note  that  since  p is  independent  of  the  pressure,  the  Gibbs  free  energy 
may  be  expressed  in  the  form 

g{T,p,<t>)  = ga{T,4’)+{lj^,  (49) 


14 


where  p0  is  a reference  pressure.  The  function  g0(T,  cf>)  is  assumed  to  have  the  form 


go{T,<i))  = 


« - ibH  (■  - B - cTI >“  (£) + 


(50) 


4a5  j \ ±m'  \J-m 

Here  e0  is  a constant  reference  energy,  the  heat  capacity  per  unit  mass  c and  latent  heat  per 
unit  mass  L are  assumed  to  be  constant,  Tm  is  the  melting  point  at  the  reference  pressure  p0? 
as  and  a are  positive  constants,  and  H(</>)  is  a double-well  potential,  which  we  will  assume  is 
given  by  H(</>)  = (f)2(l  — (f))2.  This  form  for  g0  is  consistent  with  an  internal  energy  which  is  a 
linear  function  of  temperature,  which  leads  to  the  classical  heat  equation  in  the  bulk  liquid  [25]. 

The  corresponding  expressions  for  the  entropy  and  internal  energy  are  then 


s=^{e°-rWi+iM+cln(£)’ 


(51) 


e — e0  + c(T  — Tm)  ~ r((f>)T  + 


1 1 
,4as  4a 


HW- 


Po 

pW 


(52) 


3 Examples 


In  this  section  planar  solidification  fronts  moving  with  constant  speed  V are  examined.  In  this 
case,  the  anisotropy  of  the  surface  energy  plays  no  role,  and  we  assume  T(V0)  = |V0|.  The 
solutions  are  assumed  to  depend  only  on  the  vertical  variable  z — z'  — Vt,  where  z'  is  measured 
in  the  rest  frame  of  the  solid,  and  the  one-dimensional  solutions  are  time-independent  functions 
of  2.  The  velocity  is  measured  in  the  rest  frame,  so  that  the  velocity  vanishes  in  the  absence  of 
driving  forces.  In  the  moving  frame,  the  governing  equations  (33)  - (36)  can  then  be  expressed 
as  conservation  laws  for  mass,  horizontal  and  vertical  components  of  momentum,  and  energy, 
together  with  the  phase-field  equation,  as 


0 

0 

0 


d 

dz 

d_ 

dz 

d 

dz 


[p{w-  V)], 


p(w  — V)u  — 
p(w  — V)w  + 


*>E 


2 

F 


[2 p{4>)  + ^(0)] 


(53) 

(54) 

(55) 
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° dz 


p(w  - V)  {e  + i(»2  + u,2)}  - ef(w  - V)  (g)  - fcWg 


- - CW)  + V4wg  + «/  jp  + ^4  (g 

0 = 4g-"(-v)S-pg. 


(56) 

(57) 


respectively.  Here  for  simplicity  we  consider  a two-dimensional  flow  u = + with  horizontal 

and  vertical  velocity  components  u and  w in  the  x and  2 directions,  respectively.  The  material 
properties  k,  p,  and  A are  assumed  to  be  constant  in  each  bulk  phase,  depending  only  on  <j>. 


3.1  Equilibrium  of  a Planar  Interface 

The  solution  for  a stationary  planar  interface  in  the  absence  of  convection  includes  the  conditions 
for  thermodynamic  equilibrium,  and  represents  a special  case  of  the  discussion  in  Section  2.2. 
The  sample  is  assumed  to  occupy  the  region  — oo  < 2 < 00,  with  </>(z)  — > 1 as  2 — » — 00  and 
0(2)  — >•  0 as  2 — >•  00.  The  density  and  horizontal  momentum  equations  are  satisfied  identically. 
The  energy  equation  admits  an  isothermal  solution,  and  the  remaining  equations  (55)  and  (57) 
give 

•'S  - ■ “•  |58> 

and 

p = p~-Y\Tz)’  (59) 

where  is  the  common  value  of  the  pressure  in  the  bulk  phases  where  the  gradient  of  the  phase 
field  tends  to  zero.  By  using  (58)  and  (59),  the  derivative  of  g(p,  T,  4>)  with  respect  to  2 vanishes, 
and  so  represents  a first  integral  for  the  system  as  noted  in  Section  2.2.  Upon  elimination  of 
the  pressure,  Eqn.  (58)  represents  a second  order  differential  equation  for  the  determination  of 
0;  some  numerical  examples  are  given  below. 

Equating  the  bulk  values  of  the  chemical  potential  by  setting  g(poo,T,0)  = g(p(X)1T\  1)  in 
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Eqn.  (49)  gives  an  integrated  form  of  the  Clausius-Clapeyron  relation  [59], 

{m> 

which  gives  the  dependence  of  the  phase  transition  temperature  on  the  pressure;  here  we  recall 
that  Tm  is  the  melting  point  at  the  reference  pressure  po.  Bulk  equilibrium  values  for  a one- 
dimensional planar  equilibrium  state  are  therefore  completely  determined  if,  say,  the  far-held 
value  of  the  pressure  in  the  solid  is  given,  since  the  far-held  pressure  in  the  liquid  is  the  same, 
and  the  Clausius-Clapeyron  relation  then  provides  the  value  of  the  temperature.  Bulk  values 
for  the  other  thermodynamic  variables  follow  from  the  known  values  of  the  temperature  and 
pressure. 

We  note  that  LaCombe  et  al.  [60]  have  proposed  to  take  advantage  of  the  relatively  rapid 
response  of  the  melting  point  to  pressure  changes  in  dynamic  solidihcation  studies.  For  example, 
during  dendritic  growth  the  response  of  the  tip  operating  conditions  to  pressure-induced  changes 
in  the  bulk  melting  point  can  be  examined  in  this  way.  We  also  note  that  Maruyama  et  al. 
[61]  have  examined  transitions  in  the  kinetic  growth  shapes  of  ice  J^,  from  a circular  disc  to 
a hexagonal  plate,  in  response  to  pressure-induced  alterations  of  the  melting  point  near  the 
roughening  transition. 

3.2  Density-Change  Flow 

Here  we  consider  the  steady  flow  normal  to  a moving  planar  solidification  front  that  is  gener- 
ated by  unequal  solid  and  liquid  densities;  the  horizontal  momentum  equation  is  then  satisfied 
identically.  For  simplicity  of  discussion  we  consider  motion  that  is  dominated  by  the  effects  of 
interface  attachment  kinetics  rather  than  diffusion;  that  is,  diffusion  is  sufficiently  rapid  that  the 
system  remains  in  thermal  equilibrium,  and  the  rate  of  solidification  is  determined  by  the  devi- 
ation of  the  temperature  from  the  equilibrium  melting  point.  We  therefore  ignore  the  equation 
for  energy  conservation,  and  assume  an  isothermal  system. 


(Poo  - Po) 


1 

Ps 


1 

PL 
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The  conservation  of  mass  equation  (53)  implies  that 


p(w  -V)  = J, 


(61) 


where  the  constant  Jm  represents  the  mass  flux  through  the  interface.  For  a stationary  solid  we 
have  that  Jm  = — psV , and 


The  flow  velocity  in  the  liquid  far  from  the  interface  is  then  w g3  = —V(ps/pL  ~ 1),  representing 
a contraction  flow  if  the  solid  is  more  dense  than  the  liquid. 

The  vertical  momentum  equation  (55)  gives 


(63) 


where  the  constant  p 1?  represents  the  pressure  in  the  solid  in  the  far  field.  The  difference  in  the 
far-field  pressures  in  the  solid  and  liquid  is  then  given  by 


(64) 


which  is  analogous  to  the  “vapor  recoil”  effect  in  a liquid- vapor  system;  i.e.,  a contraction  flow 
raises  the  pressure  in  the  solid  relative  to  that  in  the  liquid. 

In  addition  to  the  jump  in  the  bulk  pressure  for  a moving  interface,  there  is  a jump  in 
the  chemical  potential  that  is  determined  by  the  phase-field  equation  (57);  this  results  in  a 
generalized  version  of  the  Clausius-Clapeyron  equation  in  which  the  effect  of  pressure  on  the 
melting  point  is  altered  by  the  rate  of  solidification.  A calculation  gives 


(65) 


where  we  recall  from  Eqn.  (62)  that  w(z)  is  proportional  to  V and  vanishes  if  ps  = Pl • This 
shows  that  the  jump  in  chemical  potential  has  a term  linear  in  the  solidification  rate  that  depends 
on  the  phase-field  profile  <j>{z ),  a term  proportional  to  V2  involving  the  square  of  the  density 
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difference,  and  a viscous  dissipation  term  proportional  to  V (since  ( dw/dz )2  is  proportional  to 
V2).  By  using  Eqn.  (64)  for  the  pressure  difference  p ’<?  — this  expression  can  be  written  in 
the  form 


+^/j2^)+a(«](£)  dz ■ (66) 

For  the  special  case  of  pL  = ps  only  the  first  integral  on  the  right  hand  side  remains,  and  we 
recover  the  conventional  phase-field  description  of  kinetically-controlled  growth,  in  which  <j)(z) 
has  a hyperbolic  tangent  profile,  and  the  growth  velocity  is  proportional  to  the  product  of  1/M 
and  the  deviation  of  the  temperature  from  its  equilibrium  value.  The  general  case  requires  a 
numerical  solution  to  determine  the  phase-field  profile;  in  the  following  section  we  describe  an 
approximate  solution  which  is  valid  if  the  density  difference  between  the  liquid  and  solid  phases 
is  small. 


3.2.1  Asymptotic  Solution 

In  this  section  we  study  the  simplified  case  where  the  solid  and  liquid  densities  are  nearly  equal. 
We  define  the  density  mismatch  parameter,  (5,  by 

Ps  = Pl{  1 + 6).  (67) 

and  consider  the  solution  of  (53)-(57)  in  the  limit  <5—^0.  We  introduce  the  expansions 

(j>  = 0o  + $<f>i  + • • • , (68) 

V = y0  + (Wi  + ...,  (69) 

and  solve  Eqn.  (57)  for  0 and  V.  The  vertical  component  of  velocity  and  the  pressure  are  given 
by  (62)  and  (63),  respectively.  The  density  profile  is  given  by  Eqn.  (48),  and  the  resulting 
vertical  component  of  velocity  is  small,  with  w = O(S). 

The  0(1)  problem  for  </>0  is  given  by 


2dV<>  , dgo  n 

eF~TT  + MV0—, Pl—  = 0, 


dz 2 


dz 


d(j) 


(70) 
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where  here  and  below  go  is  the  free  energy  evaluated  at  </>0  (under  isothermal  conditions).  This 
equation  has  the  solution 


<t>o{z)  = 1 [1  - tanh (z/Q) , 


(71) 


where  the  constant  4,  which  characterizes  the  thickness  of  the  interface,  is 


4 — 2cp‘ 


2 a0(T) 
Pl 


(72) 


Here 


1 = 1 If  T \ 

ao(T)  ~ a + as  V TM ) ’ 

and  the  leading  order  interface  speed  is  related  to  the  temperature  by 


(73) 


MV0  — 6LeF^/2pLa0(T)  ^1  — j . (74) 

In  the  limit  of  vanishing  double-well  height  ratio,  a/as  ->  0,  the  velocity  of  the  interface 
depends  linearly  on  the  difference  in  temperature  from  the  melting  temperature  Tm-  When 
T < Tm  the  leading-order  interface  velocity  U0  is  positive  (i.e.  the  liquid  solidifies)  and  when 
T > Tm  the  interface  velocity  is  negative  (i.e.  the  solid  melts). 

The  0(5)  problem  for  (pi  is  given  by 


£<4 


2 d2(f) i 
"F  dz 2 


+ MV0 


d®i 

dz 


d2go  ^ 

pL  a</>2  ^ 


K, 


(75) 


where 


n = 


■MV^  - MV0[1  - r(<t>0)}d<t>° 


dz 


+ PLr{<l>o) 


dz 


dgo 

d<\) 


oo  2 

i 

CM 

o 

-e- 

Ps  -Po-  2€f 

[dz) 

r'W>o)- 


(76) 


Here  we  have  used  a density  profile  of  the  form  (48). 

The  operator  C appearing  in  equation  (75)  can  be  written  in  self-adjoint  form  by  first  mul- 
tiplying the  equation  by  exp (MV0z/e2F).  Upon  noting  that  dfo/dz  is  a solution  to  the  homoge- 
neous problem,  C(j>i  = 0,  we  deduce  the  solvability  condition 


0=1  °vz 


°°  e-^ndz, 


— CO 


dz 


(77) 
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where  77  = MV0/e2F.  Evaluating  this  integral  gives  that 

MViQi(V0)  = -(pf  - p0)Q2(V0)  + MV0Q3(Vo)  + ^^Q4(V0),  (78) 

where  Qi{V0) , Q2{V0) , Qs{V0) , Q4(V0)  are  the  following  integrals: 

/oo 

exP(7?C)  y>o(0)2dC,  (79) 

-OO 

/OO 

exp (rjQ  r'(<f>0(Q)(l>'0(Od(,  (80) 

-OO 

/OO 

exp(»?C)  [2t-(0o(C))  - l]K(C)]2dC,  (81) 

-OO 

Qi( V o)  = ^exp(r?C)4:  {[^(OlM^C))}  (82) 


Eqn.  (78)  gives  the  the  velocity  correction,  V1;  in  terms  of  V0  and  p™  — po.  Note  that  Vq  is 
related  to  the  temperature  by  Eqn.  (74)  and  hence  Eqn.  (78)  can  be  regarded  as  determining 
the  the  interface  velocity  in  terms  of  the  pressure  and  temperature.  This  relation  represents  the 
asymptotic  approximation  to  0(5)  of  the  nonequilibrium  Clausius-Clapeyron  relation  given  in 
Eqn.  (66).  From  Eqn.  (66)  we  observe  that  to  the  first  order  in  5 the  only  contribution  from  the 
nonequilibrium  terms  is  due  to  interfacial  kinetics  represented  by  the  first  integral  on  the  right 
hand  side. 

The  results  of  this  asymptotic  analysis  are  presented  belowr  in  conjunction  with  the  numerical 
results  for  general  values  of  the  solid  and  liquid  densities. 

3.2.2  Numerical  Solution 

In  this  section  we  present  numerical  solutions  of  the  one-dimensional  system  of  governing  equa- 
tions (53)-(57).  Our  aim  is  to  solve  for  the  phase- field  profile  and  the  solidification  front  velocity 
as  the  temperature  and  far-held  pressure  in  the  solid  vary,  and  compare  the  results  with  our 
small-5  asymptotic  analysis. 

We  solve  Eqn.  (57)  for  the  phase  held  with  w and  p given  in  terms  of  0 by  equations  (62) 
and  (63),  respectively.  The  numerical  procedure  involves  a hnite  difference  discretization  of  the 
governing  equation  (57)  on  the  domain  — L < z < L (where  L is  taken  to  be  sufficiently  large  so 
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as  to  not  influence  the  results  of  the  calculations).  Neumann  boundary  conditions  = 0) 

are  applied  at  z — ±L  and  we  require  that  the  integral  of  0 — 1/2  is  zero,  which  both  fixes 
the  position  of  the  interface  near  z = 0 and  provides  an  additional  equation  that  allows  the 
determination  of  the  interface  velocity  for  given  values  of  the  temperature  and  far-held  pressure 
in  the  solid.  The  discretized  system  of  nonlinear  equations  is  solved  using  Newton  iteration  by 
employing  the  software  SNSQ  [62]. 


o 

6/3 


(T-Tm)/Tm  x 105 


Figure  1.  A plot  of  the  dependence  of  the  interface  velocity  on  temper- 
ature for  three  different  far-held  pressures  in  the  solid.  The  data  used  is 
given  in  Table  1 and  corresponds  to  lead.  The  solid  curves  are  from  the 
numerical  solution  of  the  governing  equations  and  the  dashed  curves 
show  the  small  6 asymptotic  result  evaluated  from  (78)  for  S = 0.035. 
From  top  to  bottom,  the  curves  correspond  to  p — p0  = 17  bar,  0 bar, 
and  —17  bar,  respectively. 
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In  Fig.  1 the  solid  curves  display  numerical  results  for  the  variation  of  the  interface  velocity 
with  temperature  for  three  different  values  of  the  pressure  using  the  material  parameters 
for  lead  given  in  Table  1.  We  put 

p(4>)  = psr(<^)  + Mi[l  - r(<j>) j,  (83) 

where  = </>2(3</>  — 1),  and  ps  and  pl  are  the  dynamic  viscosities  of  the  bulk  solid  and 
liquid  phases,  respectively.  We  use  ps/ pl  = 1 and  = — 2p(<j>)/3  in  accordance  with  the 
Stokes  hypothesis.  In  these  calculations  ps  > Pl  and  so  from  the  equilibrium  form  of  the 
Clausius-Clapeyron  relation  (60)  we  expect  the  pressure  to  increase  with  temperature  which  is 
in  agreement  with  the  results  for  V = 0 shown  in  this  figure.  In  Fig.  1 we  also  observe  that  the 
interface  velocity  decreases  with  temperature  in  a roughly  linear  fashion  at  fixed  pressure.  The 
effect  of  changing  the  pressure  is  to  shift  the  V(T)  curves  while  preserving  their  slope.  This 
behavior  is  confirmed  by  our  expression  for  the  non-equilibrium  Clausius-Clapeyron  equation 
(66),  when  we  observe  that  the  integrals  on  the  right  hand  side  are  all  positive.  The  dashed 
curves  in  this  figure  represent  the  asymptotic  approximations  to  these  curves  that  are  obtained 
in  the  limit  of  small  density  differences. 

We  also  investigated  the  effect  of  increasing  the  viscosity  ratio,  ps/pl-  We  observe  from 
the  non-equilibrium  Clausius-Clapeyron  relation  (66)  that  the  viscous  dissipation  may  increase 
without  bound  as  the  solid  viscosity  is  increased  and  result  in  the  relationship  between  tem- 
perature and  pressure  being  significantly  disrupted.  However,  we  found  from  our  numerical 
calculations  that  for  the  material  parameters  given  in  Table  1 for  lead  with  a pressure  difference 
p — Po  — 0 that  this  is  a very  weak  effect;  the  melting  temperature  varied  by  4 % for  an  increase 
in  ps/pl  from  1 to  106.  For  a greater  interface  thickness  of  10-4  cm  the  change  in  the  melting 
temperature  is  2 x 10-3  % for  the  same  increase  in  ps/pi . 
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Table  1:  Properties  used  in  Calculations  [64>65] 


Property 

Value 

Units 

Liquid  thermal  conductivity 

Kl 

0.159 

J/(cm  s K) 

Solid  thermal  conductivity 

Ks 

0.297 

J/(cm  s K) 

Liquid  thermal  diffusivity 

kl 

0.108 

cm2/s 

Solid  thermal  diffusivity 

Ks 

0.202 

cm2/s 

Melting  point 

Tm 

600 

K 

Heat  of  fusion 

Lv 

256.0 

J/cm3 

Kinetic  coefficient 

P 

33 

cm/ (s  K) 

Kinematic  viscosity 

V 

2.43  x 10~3 

cm2/s 

Liquid  density 

PL 

10.66 

g/cm3 

Interface  thickness 

4 

1.0  x 10“8 

cm 

Density  change 

( Ps/pL ) - 1 

0.035 

- 

3.3  Shear  Flow 

We  consider  the  case  of  a shear  flow  parallel  to  a planar  interface  that  propagates  with  a constant 
velocity  V.  Far  from  the  interface,  in  the  liquid,  we  assume  that  the  component  of  fluid  velocity 
parallel  to  the  interface  is  JJ^.  We  take  a coordinate  system  coincident  with  the  moving  interface, 
which  is  given  by  the  plane  z = 0.  For  simplicity  we  consider  the  situation  where  the  system  is 
isothermal,  the  density  of  the  both  the  solid  and  liquid  phases  are  equal,  the  surface  energy  is 
isotropic,  and  the  dynamic  viscosity  has  the  form 


n(<t>)  = Hs<f>  + (1  - <t>)HL- 


(84) 


Under  these  assumptions  the  velocity  field  is  given  by  u = u(z)i  and  the  governing  equations 
(53)— (57)  reduce  to 

du 


d 

dz 

d 

dz 


pVu  + /i(</>) 


dz 


P+^i 


d(j) 

dz 


0 = 6 


2 d2<t)  dg 

?Is+MVTz-pW 


(85) 

(86) 
(87) 


where  p is  the  common  value  of  the  bulk  densities.  The  phase-field  equation  is  identical  to  the 
leading  order  equation  for  0O  in  the  previous  section.  Its  solution  is  therefore  given  by  (71)  with 
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the  requirement  that  the  interface  velocity  is  related  to  the  temperature  in  the  same  manner  as 
given  in  Eqn.  (74)  [when  V is  identified  with  V0}.  The  momentum  equations  may  be  integrated 
directly  to  give 


p{z)  = Poo  ~ 


64  a0(T) 


sech4  ( — 


and 


u(z) 

Ur*. 


— I — A exp 


1 + m\  z 


2m 


)-3{ 


cosh 


z . (\  — m 

— + tanh  

L Vl  + mj. 


where  £ — pL/{pV ) is  the  viscous  boundary  layer  thickness  in  the  liquid, 

(1  — m)  4 


P = 


2m  t ’ 


(88) 


(89) 


(90) 


and  m — ps/ Pl-  The  constant  of  integration,  A,  is  chosen  to  satisfy  the  condition  that  u(—Z)  = 
0,  where  Z > 0: 


A = exp 


/I  + m\  Z' 
V 2m  ) £ . 


+ tanh  1 


1 — m 
1 + m 


(91) 


Taking  the  sharp  interface  limit,  £e/£  — ► 0,  it  is  found  that  u(z)  ~ U00{l—exp[—(z+Z/m)/£]}, 
for  z > 0 (in  the  liquid),  and  u(z)  ~ bfoojl  — exp[— (z  + Z)/(£m)]},  for  z < 0 (in  the  solid). 
Subsequently,  taking  the  limit  in  which  the  viscosity  of  the  solid  is  much  greater  than  that  of 
the  liquid,  m — >•  oo,  it  is  found  to  leading  order  in  the  liquid  that  u(z)  ~ Uoo{l  — exp (—z/£)}, 
and  in  the  solid  that  u(z)  ~ U00(z  + Z)/(£m ),  i.e.,  a linear  shear. 

Thus  in  the  sharp  interface  limit  with  the  viscosity  of  the  solid  much  greater  than  that  of 
the  liquid  we  recover,  at  leading  order,  the  exact  solution  for  u(z)  in  the  liquid  to  the  underlying 
free  boundary  problem,  i.e.,  the  asymptotic  boundary-layer  profile  u(z)  = Uoo{l  — exp (—z/£)}. 
In  the  solid  phase  we  obtain  a uniform  shear  flow,  with  a strain  rate  of  magnitude  Uool{2£m). 
The  velocity  at  the  center  of  the  interface  (z  = 0)  is  U^Z /(£m)  and  so  we  find  that  the  no-slip 
condition  at  the  interface  is  almost  satisfied,  the  magnitude  of  the  error  decreasing  in  a manner 
inversely  proportional  to  m. 

Plots  of  u(z)  for  different  values  of  m with  £J£  — 0.01  are  given  in  Fig.  2.  For  all  values 
of  m we  observe  a thin  interfacial  layer  separating  the  linear  shear  in  the  solid  from  the  shear 
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flow  in  the  liquid.  We  see  that  as  m increases  the  flow  in  the  solid  region  is  suppressed  and 
adopts  a decreasingly  small  uniform  shear  while  the  flow  in  the  liquid  approaches  the  asymptotic 
boundary-layer  profile,  consonant  with  our  asymptotic  analysis.  The  no-slip  condition  at  the 
interface,  z = 0,  is  satisfied  with  greater  accuracy  as  m increases. 


o 


ing  equations  for  le/i  = 0.01  and  three  different  values  of  the  viscosity 
ratio  hs/hl-  From  top  to  bottom,  the  curves  correspond  to  Hs/hl 
equal  to  10,  30  and  500. 

The  above  analysis  results  in  a closed-form  solution  for  the  specific  choice  (84)  of  the  dynamic 
viscosity  //(</>);  other  choices  are,  of  course,  possible,  and  Eqn.  (84)  was  chosen  as  a compromise 
that  maintains  consistency  with  the  form  (83)  used  in  the  previous  section  while  allowing  a 
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simple  exact  solution.  In  general,  the  detailed  behavior  of  the  fluid  flow  in  the  interfacial  region 
would  be  expected  to  depend  on  the  specific  form  chosen  for  /i(0).  For  example,  interpolating 
by  the  reciprocal,  as  in  1 / ii(<j>)  — </>/ ns  + {1  — (f>)/nL,  also  results  in  a closed-form  solution. 
Similarly  the  viscous  dissipation  in  the  interface  will  also  depend  on  the  form  of  //(</>).  This 
in  turn  will  influence  the  degree  to  which  the  melting  temperature  is  influenced  by  the  ratio 
Us/ Hl  discussed  above  for  a density-change  flow.  We  note  that  in  a related  situation  in  which  a 
diffuse-interface  model  was  used  to  study  non-equilibrium  effects  during  directional  solidification 
of  a binary  alloy  [66],  the  details  of  some  interfacial  quantities,  such  as  a characteristic  trapping 
velocity,  were  found  to  be  sensitive  to  the  specific  form  chosen  for  interpolating  the  diffusivity 
through  the  interface;  such  considerations  likely  apply  to  the  present  case  as  well. 

4 Conclusions 

In  this  paper  we  have  developed  a phase-field  model  for  the  solidification  of  a pure  material  that 
includes  convection.  The  model  has  been  derived  using  the  formalism  of  irreversible  thermody- 
namics, employing  gradient  energy  and  entropy  terms  that  incorporate  the  effects  of  capillarity 
within  the  framework  of  a diffuse  interface  model.  Our  solidification  model  has  two  distinctive 
aspects:  we  treat  both  the  solid  and  liquid  phases  as  Newtonian  viscous  fluids  (with  a high 
ratio  of  solid  to  liquid  viscosity,  and  the  phase  boundaries  are  endowed  with  anisotropic 

surface  energy  by  using  the  generalized  ^-vector  formulation.  Both  features  are  consonant  with 
our  intention  to  model  a solid-liquid  system.  Our  formulation  neglects  elastic  effects  in  the 
solid.  We  work  with  pressure  and  temperature  as  independent  thermodynamic  variables  which 
permit  a quasi-incompressible  formulation  in  which  the  densities  can  be  uniform  in  each  phase. 
This  allows  us  to  model  volume-change  effects  on  solidification  which  we  illustrated  for  the  case 
of  constant  velocity  isothermal  unidirectional  solidification.  In  this  setting  wTe  derived  a non- 
equilibrium form  of  the  Clausius-Clapeyron  relation  which  describes  how  the  relation  between 
pressure  and  temperature  is  modified  by  the  motion  of  the  interface.  The  ability  of  this  two- 
fluid  model  to  approximate  the  no-slip  condition  at  the  interface  was  examined  by  considering 
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a model  shear  flow  parallel  to  a planar  interface.  In  this  situation  we  found  an  exact  solution 
which  showed  that  the  correction  to  the  no-slip  boundary  condition  is  first  order  in  the  quantity 
Hl/^s  f°r  a simple  model  in  which  the  viscosity  depends  linearly  on  the  order  parameter  (j>. 

We  have  focused  our  attention  on  the  derivation  of  the  model  and  some  important  test  cases 
involving  planar  interfaces.  In  future  work  we  plan  to  investigate  the  sharp  interface  limit  for 
curved  interfaces  as  well  as  use  the  model  as  a computational  vehicle  for  the  investigation  of 
flow  effects  in  dendritic  solidification,  including  the  effects  of  buoyancy  forces. 
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Appendix 

We  have 


If 

dt  Jn(t) 


\vvdv 


'n{t) 

f 

'n(t) 

f 

'n{t) 


iar2  i ,_2_ 

+ -v  ■ r 2u) 

2 dt  2 v ' 


dV, 


r^  + £-(rvr)  + ir2v-£ 

dt  v ’ 2 


dV, 


dV, 


(92) 


where  we  have  used  the  fundamental  relation  dT(p)  = £ • dp,  with  p = V</>.  We  note  that 


v'(r*75t)  = ^v-(r|)  + rf  + rvc  = * ® v*.  (93) 

where  the  tensor  A = £ ® V<j>  has  components  Ajk  = £jd(/)/dxk,  and  the  double  contraction  of 
the  tensor  product  is  denoted  by  Vu  : A = dukjdx3  A 3k.  This  gives  that 


-/ 

dt  Jn(t) 


\r2(v<p)dv 


'n(t) 
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_ Dd 
Dt ) Dt 
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